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Abstract 

Tumor growth has long been a target of investigation within the context of math- 
ematical and computer modelling. The objective of this study is to propose and 
analyze a two-dimensional probabilistic cellular automata model to describe avas- 
cular solid tumor growth, taking into account both the competition between cancer 
cells and normal cells for nutrients and/or space and a time-dependent proliferation 
of cancer cells. Gompertzian growth, characteristic of some tumors, is described and 
some of the features of the time-spatial pattern of solid tumors, such as compact 
morphology with irregular borders, are captured. The parameter space is studied in 
order to analyze the occurrence of necrosis and the response to therapy. Our findings 
suggest that transitions exist between necrotic and non-necrotic phases (no-therapy 
cases), and between the states of cure and non-cure (therapy cases). To analyze 
cure, the control and order parameters are, respectively, the highest probability of 
cancer cell proliferation and the probability of the therapeutic effect on cancer cells. 
With respect to patterns, it is possible to observe the inner necrotic core and the 
effect of the therapy destroying the tumor from its outer borders inwards. 
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1 INTRODUCTION 

Neoplastic diseases are the cause of 7 million deaths annually or 12of deaths 
worldwide pp. Mathematical and computer modelling may lead to greater 
understanding of the dynamics of cancer progression in the patient [2] [3]. 
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These techniques may also be useful in selecting better therapeutic strategies 
by subjecting available options to computer testing (in silico). Continuous 
models have been proposed to describe the stages of tumor growth since the 
middle of the 20th century [3]. Initially, a tumor grows exponentially (linear 
rate). After this transient stage, the growth rate decreases and a steady state 
is attained, due to several factors including a lack of nutrients and hypoxia. 
This nonlinear behavior characterizes avascular tumor growth when neovascu- 
larization has not yet been triggered. The decelerating avascular growth may 
be guided by different rules such as, for example, Gompertzian and logistic 
functions. Gompertzian growth has been one of the most studied decelerat- 
ing tumor growth over the past 60 years [5], [6] [7]. It is found, for example, 
in some solid tumors such as breast carcinomas [8] [9]. It is also observed in 
tumors in vitro [TU] [TTj . 

Although continuous models are capable of describing the behavior of tumor 
growth, it would appear more reasonable to adopt a discrete approach when 
describing the prevascular stage. Due to the fact that the angiogenic process 
has not triggered early tumor growth, few cancer cells are present and growth 
depends predominantly on the interactions of these cells with adjacent cells 
and with the environment [12]. In addition, in the discrete approach, it is 
easier to capture the time-spatial pattern generated by the model in order to 
compare it with actual patterns [13] . Some cellular automata models [H] [15] 
and hybrid cellular automata [TJ] [IS] [12] have been proposed to study tumor 
growth. 

Another important topic that is analyzed in mathematical models is the re- 
sponse to therapy, including how tumor growth changes under the effect of 
a drug [18]. The focus is directed towards identifying the optimal therapy to 
maximize the effect on cancer cells and minimize the effect on normal cells 
[19] [20] . Although various continuous chemotherapy models exist [21], [22], 
[23] . to the best of our knowledge the majority of the discrete models cited in 
the literature have not yet been used to investigate this topic. 

The objective of this work is to propose a two-dimensional cellular automata 
model consisting of 4 states (empty site, normal cell, cancer cell or necrotic 
tumor cell) to describe avascular tumor growth. It should be emphasized that 
in this study the term tumor growth is used to refer to the number of cancer 
cells rather than the volume of the tumor; in other words, it is assumed that 
the tumor volume is proportional to the number of cancer cells [3]. Assuming 
that the angiogenic process has not yet been triggered, there is no increase 
in nutrients, which are uniformly distributed over the lattice. In this simple 
model, some relevant processes involved in the prevascular phase of tumor 
growth are assessed: a dynamic proliferation of cancer cells and the competi- 
tion between normal cells and cancer cells for nutrients and/or space. Since 
necrosis is often present in the prevascular stage of tumor growth [TO] [23], 
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the possibility of necrosis in the model must also be taken into consideration. 
Finally, the effect of therapy is included in order to investigate whether the 
system evolves to a state of cure. 

This paper is organized as follows. In section [2j the model is presented, to- 
gether with its local rules, parameters and the scope of the algorithm. Section 
[3] describes the simulated time series of cell density in the presence or absence 
of treatment, and shows the features of the time-spatial pattern of simulated 
solid tumors. In section [4], the parameter related to the process of necrosis and 
the effect of therapy is analyzed. Finally, in section |5j our results are discussed 
from the point of view of the phenomenon and some concluding remarks are 
made. 



2 THE MODEL 



We propose a two-dimensional (L x L) cellular automata model [25] under 
periodic boundary conditions, using a Moore neighborhood with a radius of 1. 
At the initial condition (t — 0), there is only one cancer cell (to ensure better 
visualization, this was taken from the center of the lattice). Since the intention 
it to model a non-viral tumor, normal cells would not be transformed into 
cancer cells with the exception of the cancer cell that triggers tumor growth 
at t = [26J. The lattice represents a tissue sample; there is a cell in each site 
that may be in one of four states: normal cell (NoC), cancer cell (CC), necrotic 
cell (NeC) or empty site (ES). We assume that the nutrients are uniformly 
available over the lattice. In this respect, lack of space is identified with lack 
of nutrients in our model. 

There is a growth potential P c value associated with each normal or cancer 
cell. Although the two-dimensional character of the model mimics the in vitro 
situation, the growth potential of cancer cells simulates the three-dimensional 
tumor in vivo in the sense that it represents the total number of cancer cells, 
i.e., in addition to the cancer cells on the lattice, the cancer cells generated by 
these cells. 



The local rules are such that: 



(i) The initial value of the mitotic probability of cancer cells is represented 
by a parameter p that measures the available resources at the beginning 
of the tumor. After that, it decreases by a factor Ap mito t until reaching the 
null value: 



Ap. 



mitot 



exp 



4 n cc (t) , 



(1) 
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in which n noc and the number of normal and cancer cells at time 

t, respectively. As shown, Ap mitot (t) depends only on the dynamics, setting 
up a feedback inhibition mechanism [T3]: as the tumor grows, Ap mitot (t) 
decreases because of the combined effect of the decrease in the number of 
normal cells and the increase in the number of cancer cells. In order to in- 
tensify the effect of this mechanism (see [27]), an exponent 2 in equation 
([I]) is considered. Since there is no new available source of nutrients and/or 
space, it decreases as the density of cancer cells increases because the avail- 
able nutrients and/or space are reduced. The effect of the proliferation of 
cancer cells is that their growth potential increases by a unit at each time 
step. 

(ii) The cancer cells compete with normal cells for the empty sites, depend- 
ing on the potential growth of neighboring cells. According to the majority 
rule, a normal cell is displaced by a cancer cell following local battles oc- 
curring between healthy and cancerous cells [15]. 

(iii) If the growth potential of a cancer cell reaches a threshold value that is 
a fraction / of the lattice size L, it becomes necrotic and its growth potential 
falls to zero. 

(iv) Both normal and cancer cells may die, with probabilities Pdrugn and 
Pdmgc, respectively, due to the continuous infusion of a drug that is applied 
after tap time steps; in this case, the site becomes empty. 

(v) if there are no cancer cells in the neighborhood of a dead cell (empty 
site), regeneration of normal cells occurs; if the cancer cells in the neigh- 
borhood of a necrotic cell die as a result of the therapy, the necrotic cell is 
eliminated. 

The algorithm was computationally implemented in FORTRAN 77 in accor- 
dance with the following steps: input data; calculate Ap mito t(t)] identify the 
state of the cell (choose one of the subroutines: normal cell (NoC), cancer 
cell(CC), necrotic cell(NeC), empty site (ES); update the cells of the lattice; 
after N iterations, output data. 

The input data are the following cellular automata (CA) parameters: 

1) The spatial parameters: lattice size L; necrosis threshold fraction / of 
lattice size; 

2) The temporal parameters: the length of the time series t final and the 
initial time of therapy infusion t ap ; 

3) The probabilities of: the initial proliferation of cancer cells po~, the effect 
of therapy on normal cells and cancer cells {pdrugn and Pdrugc)- 

The output data are the time series of the density of each type of cell and the 
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final configuration of the lattice at any time step. In addition, the time-spatial 
configurations, controlled by a package denominated g2 [2H] whose commands 
are inserted into the computer program in FORTRAN, are generated in "real 
time". This package may be used in C, PYTHON and PERL. 

The following is a description of each subroutine based on the local rules: 

a) Normal Cell (NoC) - a random number y is compared to Pdrugn- If 
y < Pdrugn, the growth potential P noc is confirmed: if P noc = 0, the cell dies 
and the site becomes empty; otherwise, it remains occupied by a normal cell 
but P noc — 0. If y > Pdrugn and if there is at least one neighboring cancer 
cell, then the normal cell is 'dislocated', P noc = and the site becomes 
empty; otherwise it remains a normal cell. 

b) Cancer Cell (CC) - if all of its neighbors are cancer cells and its growth 
potential reaches a fraction / of lattice size L, the cancer cell becomes 
necrotic. Otherwise, a number y is randomly chosen. If y < pdmgc and 
t > t ap , the potential P cc is confirmed: if P cc > 0, it is reduced by a unit and 
the cell remains a cancer cell; otherwise the cell dies and the site becomes 
empty. Finally, if y > Pdrugc, the cell remains a cancer cell; however, its 
growth potential P cc increases by a unit. 

c) Necrotic Cell (NeC) - if at least one of its neighbors is neither a 
cancer cell nor a necrotic cell, it is eliminated and the site becomes empty; 
otherwise, it continues necrotic. 

d) Empty Site (ES) - if there are cancer and normal cells in its neigh- 
borhood, the local battle between cancer cells and normal cells is such that 
if the sum of the potential growth of its neighboring cancer cells is greater 
or equal to the sum of the potential growth of its normal cell neighbors, 
then the empty site is occupied by a cancer cell that diffuses from one of 
the randomly chosen neighbors; otherwise, it is occupied by one of the ran- 
domly chosen normal cells that were previously dislocated. If there are only 
normal cells in its neighborhood, it becomes a normal cell through a process 
of regeneration. Finally, if none of its neighbors are cancer cells or normal 
cells, it remains empty. 



3 RESULTS: SIMULATED TIME SERIES AND TIME-SPATIAL 
PATTERNS 



Computational simulations of the model were performed in order to analyze 
two classes of behavior: cases in which no treatment was given and treatment 
cases. In each subsection, the time series of the simulated tumor as well as 
time-spatial patterns are shown. 
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Time Step Time Step 

Fig. 1. Nonnecrotic tumor: the average of simulated time series of the density 
of: a) normal cells, b) cancer cells. We consider M samp i es = 200 and the following 
parameter values: L = 251, po = 0.95, / = 0.6, Pdrugn = Pdrugc = 0.0. 



a) b) 



180000- 




ti me step Time Step 

Fig. 2. Non-necrotic case: assuming M samp i es = 200 and the same parameters values 
of Figure the average of simulated time series of: a) the number of cancer cells 
(grey color); b) the growth potential of cancer cells (black color). The Gompertzian 
fitting is applied to (a) (black color) those in which parameters of time series are 
a = (6.58 ± 0.09) x 10~ 2 , (3 = (1.631 ± 0.008) x 10" 2 and n = (1.6 ± 0.06) x 10 2 . 



3.1 NO TREATMENT CASE 



In the case of no treatment, the following parameters are considered: Pdrugn — 
Pdrugc = tap = 0. For fixed values of L, t final and po, but different values of 
/, in Figures [T] and |3j the time series of the density of cells for nonnecrotic 
and necrotic tumors, respectively, are shown. In both cases, the cell densities 
reach saturated values due to the effects of the competition between normal 
cells and cancer cells, and the time-dependent mitotic probability. Comparing 
Figures [T]d and[3]3, the stationary value of cancer cell density is clearly greater 
in necrotic tumors than in nonnecrotic ones. This is a consequence of the 
fact that Apmitotit) assumes smaller values in necrotic tumors compared to 
nonnecrotic ones because it does not depend on the density of necrotic cells. 
In both cases, the average of the different samples is considered, corresponding 
to different seeds of random numbers. 



6 




1 DO 200 300 400 500 600 700 

Time Step 




!0a 200 300 400 500 600 700 

Time Step 



Fig. 3. Necrotic case: the average of simulated time series of density of a) normal 
cells, b) cancer cells. We consider M samp i es = 200 and the following parameter 

values: L = 251, p = 0.95, / = 0.2, Pdrugn = 0.0, Pdrugc = 0.0. 
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Fig. 4. Assuming M samp i es = 200 and the same parameters values of Figure 
[3j the average of simulated time series of: a) the number of cancer cells (grey 
color); b) the growth potential of cancer cells (black color). The Gompertzian fit- 
ting is applied on (a) (black color) those in which parameters of time series are 
a = (5.86±0.06) x 10~ 2 , (3 = (1.450±0.006) x 10~ 2 and n = (1.926±0.005) x 10 2 . 

Figures [2^i and |4^i show that the tumor growth obeys the Gompertzian func- 
tion both in nonnecrotic and necrotic cases. This behavior is observed with 
respect to the number of cancer cells for a range of values of po- 111 Gompertzian 
growth, the specific growth rate of the number of cancer cells decreases loga- 
rithmically: 



1 dn cc n , ( n cc 

-jr = "o - P log — 

n cc at \ n 



(2) 



where 

a) no is the initial population of cancer cells; 

b) a is the specific growth rate of n cells at t = 0; 

c) j3 measures how rapidly the curve departs from a singular exponential 
and curves over, assuming its characteristic shape. 
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The solution of ^ is 

riccif) = n exp <^[1- exp(-(3t)] \ . (3) 



It is evident that the stationary value of n cc is n CCoo = n exp(a>o//3). Accord- 
ing to the Gompertzian fitting represented by equation [3j the results of the 
simulations shown in Figures [2] and [4] correspond respectively to the following 
parameters: 

I) The nonnecrotic tumor: a = (6.58 ± 0.09) x 10~ 2 , (3 = (1.631 ± 0.008) x 
lO- 2 and n = (1.60 ± 0.06) x 10 2 

II) The necrotic tumor: a = (5.86±0.06) x KT 2 , (3 = (1.450±0.006) x 10~ 2 
and n = (1.926 ± 0.005) x 10 2 

Comparing the above parameters of (I) and (II), the behavior of nonnecrotic 
and necrotic tumors was found to be very similar. It was also found that the 
number of cancer cells obeys the Gompertzian fitting. 

An important confirmation of our model is shown by comparing the Gom- 
pertzian fitting parameters of simulated tumors with the corresponding pa- 
rameters of actual tumors [6] [?]. For instance, in the case of the testicular 
tumors shown in Table 3 of reference [6], the f3 values are in the range of 
[0.005; 0.016] day" 1 . Our simulated (3 values are within the above range if we 
consider the time step of our simulations to be one day. 

In relation to the parameters «o and no, the simulated values are not compa- 
rable to the actual values, since no information on the initial size of the tumor 
was included in our model. Both a and n are strongly dependent on that 
information. 

Finally, it would be very interesting to discover whether it is possible to re- 
late the cellular automata no-therapy parameters L, f and po with the Gom- 
pertzian parameters for the range of values of po and to assess whether the 
number of cancer cells obeys a Gompertzian growth pattern. In this case, some 
preliminary conclusions may be drawn with respect to n^: it decreases in ac- 
cordance with the necrotic parameter / but it increases linearly as a function 
of lattice size L and, exponentially with pq. Concerning the other Gompertzian 
parameters, the answer to this question is not so simple. Following these con- 
clusions leads us to a much more interesting quandary if we want to make 
the model more realistic: to estimate realistic ranges of the CA no-therapy 
parameters L, f and po based on the Gompertzian parameters of actual solid 
tumors. 

With respect to the time-spatial patterns of the simulated tumors, Figure [5] 
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Fig. 5. The spatial distribution of states at consequent time steps (see the time 
step bar) using the same parameters values of Figure [3] except L = 81 and p$ = 0.8. 
The final time step is tfi na i = 600. Dark colors (light grey, grey, and dark grey) 
correspond, respectively, to normal, cancer, and necrotic cells. 

shows the lattice configuration at some time during the steady state in the 
necrotic (Figure^) cases. Using the g2 package, it was found that, although 
tumor growth leads to the compact shape that is characteristic of solid tumors, 
the growth process is such that its boundary is irregular for any time t, as, for 
example, for the time steps represented in Figure [5j A similar type of behavior 
is observed for the time-spatial patterns in the nonnecrotic tumor. In relation 
to the process of necrosis, the necrotic region was found to be inside the tumor 
at any time t [2S] P2], as shown in Figure [5] 

3.2 TREATMENT CASE 

Figures [6^, and [6]d show the simulated time series of the number of cancer 
cells and normal cells for the different values of Pdrugc > Po an d Pdrugc < Po 
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Fig. 6. a) Successful treatment: the average of simulated time series of the number of 
normal cells (-) and cancer cells (circle) assuming M samp i es = 200 and the following 
parameters values: L = 251, po = 0.8, / = 0.2, Pdrugc = 0-9, t ap = 500, and two 
values for Pdrugn- 10 -4 (dark grey), and 10~ 2 (light grey); b) Unsuccessful treatment: 
the time series of the number of normal cells (-) and cancer cells (circle) assuming 
M samp i es = 200 and the following parameters values: L = 251, po = 0.8, / = 0.2, 
Pdrugc = 0.7, t ap = 500, and two values for Pdrugn- 10~ 4 (dark grey), and 10 -2 (light 
grey). 



corresponding to successful treatment (cure) and unsuccessful treatment (non- 
cure), respectively. In each figure, two values of Vdrugn are considered. It can be 
clearly seen that, for both values of Pdrugn, the success of the treatment does 
not change. It is also clear that the tumor size increases until tap (see figures 
[6p, and^). However when Pdrugc > Po, the reduction in tumor size starts, as 
expected, after t ap time steps. 



When treatment is successful (pdrugc > Po), it was found that, for an upper 
value of Pdrugn (Pdrugn = 10~ 2 ; light grey), more normal cells are eliminated 
compared with a lower value of Pdrugn - Pdrugn = 10~ 4 as shown in Figure 
[6^, (dark grey). Therefore, we may conclude that very small values of Pdrugn 
correspond to optimal therapy. However, in the case of unsuccessful treatment 
(pdrugc < Po), when Pdrugn is increased, both a decrease in the number of 
normal cells and a slower rate of increase of cancer cells is found (see fig (6)3). 



A systematic analysis of the parameter space presented in the next section 
will provide a more precise conclusion about the role of Pdrugn parameter. 



In order to analyze the effect of therapy on cancer cells, the time-spatial dis- 
tribution of the simulated tumors was followed using the g2 package. Figure [7] 
shows the configuration of the lattice at different time steps from the beginning 
of therapy until the tumor is eliminated. 
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Fig. 7. The spatial distribution of states at consequent time steps (see the time step 
bar) using the same parameters values of Figure [6^, except L = 81. The final time 
step is t fi na i = 1400. Dark colors (light grey, grey, dark grey, and black) correspond, 
respectively, to normal, cancer, necrotic cells, and empty sites. 



4 RESULTS: PARAMETER SPACE 



Analysis of the parameter space is important in order to confirm the robustness 
of the model and to identify the most relevant parameters for the dynamics 
of the model. 

According to the relevance of some CA parameters to the features of the 
phenomenon, this analysis of parameter space was divided into two parts: 
the occurrence of necrosis (no treatment cases) and reaching a state of cure 
(treatment cases). 
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4.1 THE OCCURRENCE OF NECROSIS 



In the first part, the parameters are again established as: Pdrugn = Pdrugc = 
tap = with the aim of evaluating the effect of necrosis, and the minimum 
value of f is obtained for each pair of values (L;p ). The values of L are 
presumed to range from 101 to 251, increasing the interval by AL = 25. 
With respect to po, the whole interval from 0.1 to 0.9 is taken into account, 
increasing Ap mitot (t) = 0.1. Figure [8] shows that a transition exists between 
the nonnecrotic (lower) and necrotic (upper) regions of parameter space. 

It was found that: 

f min = a(L)p + b(L) (4) 

where the linear and the angular coefficients, b(L) and a(L) are different for 
distinct values of the lattice size. 




Fig. 8. a) Parameter space in the case of no treatment: f mm x po X L; the following 
parameters values are fixed: Pdrugn = 0, Pdrugc = 0; tfi na i = 5000. b) the dependence 
of the angular coefficient of f m i n , a (L), with L. 

The angular coefficient a(L) = clqL 1 where a = 1.12 and 7 = 0.29 (see figure 




Since the lattice size is an intrinsic parameter of the model, that transition 
is such that the the maximal growth probability po is the control parameter, 
while the necrotic parameter / is the order parameter. 

4.2 REACHING THE CURE STATE 

If the tumor is submitted to systemic treatment represented by a probability 
Pdrugc 7^ starting at t ap ^ that affects normal cells with a probability 
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Fig. 9. Parameter subspace of treatment case: pdrugc x Po x tap', the following pa- 
rameters values are fixed: L = 251, / = 0.6, Pdrugn = 10~ 4 ; tf ina i = 5000. 

Pdrugn < Pdrugc, our aim is to establish the minimum value of Pdrugc in order to 
achieve a state of cure, i.e., no cancer cells. 



The analysis is now more complex than the one performed in subsection |4.1| in 
the sense that there are 4 important parameters that control the behavior of 
the therapeutic effect: Pdrugc, Pdrugn, tap, and p . The relevance of p is evident 



in subsection 3.2 This analysis is divided into two parts. 



In the first part, motivated by the behavior observed in figure |6^, the pa- 
rameter is defined as Pdrugn, assuming very small values (pdrugn = 10~ 4 ) to 
simulate optimal therapy. Lattice size is also established as L = 251 and the 
necrotic parameter / = 0.6. Analogously to the method applied in subsection 
4.1 , the minimum value of Pdrugc was obtained for each pair of values (p ; t ap ). 
The range of values for p and tap are, respectively, [0.1; 0.9] and [0;tf ina i\. 
Figure [9] shows the relevant role of po dividing the parameter subspace into 
two regions: the lower is the non-cure state while the upper corresponds to 
the state of cure. In this case p is again the control parameter but Pdrugc is 
the order parameter. 

In the second part, a value of p$ = 0.8 remains fixed, while the pair varies 
(Pdrugn, t ap ), again in order to obtain the minimum value of Pdrugc that would 
be sufficient to effectively eliminate the tumor, now related to the effect of the 
drug on the normal cells. The whole interval of Pdrugn was taken into account 
from very small values 10~ 4 until 0.9. It would be necessary to extend tf ina i 
in order to maintain the duration of application, since in this analysis the 
parameter t ap varied. Since this behavior is similar for any value of Pdrugn, in 



figure 10 the minimum value of Pdrugc is shown to increase as a function of t ap 



in accordance with a Lorentzian function: 

2A 1 A 3 



ram _ A , 

Pdru9C ° * A{t ap -A 2 y + Al 
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Fig. 10. Parameter subspace of treatment case: p™ugc x *»pi 
values are fixed: L = 251, / = 0.6, p = 0.6; pdrugn = 10~ 4 ; 5 t = tf ina i 
The Lorentzian fitting (grey color) parameters are Aq = 0.60, A\ = 186.61, 
A 2 = 4829.92, A 3 = 776.29. 




Fig. 11. Analysis of cure time in the treatment case: t cure xp™ugc x * a P> * ne following 
parameters values are fixed: L = 251, / = 0.6, po = 0.6, Pdrugn = 10~ 4 , tfi na i = 500. 



This result means that the minimal rate of infusion of the drug to eliminate 
the tumor has to be greater if the treatment begins later. It emphasizes how 
important it is to initiate treatment as early as possible in order to reduce the 
infusion rate of the drug. 



Finally in the cases in which the tumor is eliminated (pdrugc > Po), for fixed 
values of L, f, p and Pdrugn, the cure time (t cure ) is estimated in relation 



to the initial time application t ap , and Pdmgc (see Figure 11). Note that t cure 
is not a parameter but a consequence of the time evolution of the system. 
Figure refFigurell shows that t cure increases linearly with t ap with an angular 



coefficient equal to 1. Figure 11 also shows that, for a fixed t ap , t cure decreases 



with Pdmgc- This result also shows how important it is to start treatment as 
early as possible to reduce the amount of time required to reach the state of 



cure. 
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5 DISCUSSION AND CONCLUDING REMARKS 



The model proposed in this study is capable of capturing the Gompertzian 
behavior of avascular tumor growth. The competition between normal and 
cancer cells and the dynamic character of the mitotic probability are the rel- 
evant components of the success of this model. 

The number of cancer cells simulates tumors in vitro due to the two-dimensional 
character of the model, and their potential growth simulates a tumor in vivo 
due to its three-dimensional nature. Figures [2] and [4] show that the model is 
able to capture the dynamics of both in vivo and in vitro avascular tumor 
growth. The simulated values for the most important Gompertzian parame- 
ter j3, which characterizes the Gompertzian shape, are compatible with the 
parameter values of some tumors [6] [7j. 

The model is also able to capture necrotic and nonnecrotic tumors depending 
on the values of the parameter /. It is well-known that necrosis, unlike apop- 
tosis, is a typical phenomenon found in a group of cells that is simulated in 
our model by the changing of the potential growth of the cells. 

The time-spatial patterns reveal a tumor with a compact shape and irregular 
boundaries, as occurs in some solid tumors [12] , [30] . Evidence was also found 
to confirm the three stages of avascular tumor growth [21]: 

a) Stage I - when the tumor grows exponentially due to available resources 
(nutrients and oxygen) - see first stage in Figure [5] 

b) Stage II - when the stabilization of Ap mito t starts but there are still 
enough resources to ensure that necrosis does not occur - see second stage 
in Figure [5]d; 

c) Stage III - when necrosis may occur depending on the value of the pa- 
rameter / because the resources are insufficient to provide for tumor growth 
- see third stage in Figure in the case of necrosis. 

The next stage, the angiogenic phase of the tumor [31J, which has not yet 
been dealt with in our model, corresponds to vascular growth. To be able to 
evaluate this stage, the process of angiogenesis would have to be taken into 
account [24] . 

In the case of no treatment {pdrugc = Pdrugn — 0), the minimum value of / for 
necrosis to occur is governed by the equation Q that was obtained from the 
investigation of the parameter space. This means that, for a simulated tissue 
of dimension L, the occurrence of necrosis depends linearly on the maximum 
value of mitotic probability of the specific tumor. 

In relation to the case in which therapy was implemented, a continuous strat- 
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egy of systemic therapy, i.e. chemotherapy, was selected. Although the schedule 
of chemotherapy is usually periodic, in this case it was decided to simulate 
this less realistic situation so that the effect of the parameters p and Pdrugc 
could be compared in a simple fashion. It was thus found that, with respect 
to avascular tumor growth, when the values of the parameter pdrugn are very 
small and when pdrugc is slightly greater than po, the tumor is completely elim- 
inated (see Figure [7]). Since neovascularization has not yet been triggered, a 
state of cure may be expected to occur with a periodic schedule. 



The time-spatial patterns of the cases in which therapy was implemented (see 
Figure [8]) show that the drug acts from the borders of the simulated solid 
tumor inwards, as would be expected in the case of solid tumors. 



With respect to the effect of parameters on the state of cure, for fixed values 
of lattice size L and necrotic parameter /, and for a range of values of t ap 
and Pdrugn, the minimum value of Pdrugc again coincides with p (see Figure 
[9]). This reinforces the importance of the higher value of po in the response of 
the tumor to therapy. It retains the memory of cancer cells with respect to 
the onset of mitosis. 



Finally, Figures 10 and 11 illustrate a very relevant finding from the phe- 
nomenological point of view: the importance of initiating therapy as early as 
possible in order to reach the state of cure. The cure time was found to be 
proportional to t ap that measures the instant when the infusion starts (see 
Figure 11) and the minimum value of the therapeutic infusion to eliminate 
the tumor increases nonlinearly as a function of the starting point of therapy. 



Future studies should be carried out to generalize the model with the objec- 
tive of including the angiogenic process and the periodic schedule of systemic 
therapy. However, the most important perspective of this line of investigation 
is to compare the model with the in vitro tumor growth of cells from specific 
tissue samples and to compare parameter values. It would then be possible to 
relate the Gompertzian fitting parameters with the parameters of the model. 
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